/*
 *  mecsystem.cpp
 *  EPPI-0.0
 *
 *  Created by Juan Carlos Garcia Orden on 15/11/07.
 *
 */

#include "mecsystem.h"
#include "simulation.h"

MecSystem::MecSystem()
:	name(""),
	type(""),
	dim(3),
	mBodies(),
	mJoints(),
	mPoints(),
	childSysNamesG(),
	mChildSys(),
	parentSysName(""),
	parentProc(0),
	isWellDefined(false),
	mChangedPoints(),
	pointSets(),
	mPointIniPos(),
	mPointIniVel()
{}

// ***

MecSystem::MecSystem(string nam, ifstream &in, int di, Process *pProc, ofstream &ou )
:	name(nam),
	type(""),
	dim(di),
	mBodies(),
	mJoints(),
	mPoints(),
	childSysNamesG(),
	mChildSys(),
	parentSysName(""),
	parentProc(pProc),
	isWellDefined(false),
	mChangedPoints(),
	pointSets(),
	mPointIniPos(),
	mPointIniVel()
{
	
	char keyword[20];

	isWellDefined = false;
	
	if (parentProc->num==rankk)
	{
		ou.open(parentProc->outFileName,ios::app);
 	}

	while (in >> keyword && strcmp(keyword,"ENDSYSTEM") )
	{
		if(!strcmp(keyword,"TYPE"))
		{
			in >> type;
			if (parentProc->num==rankk)
			{
				cout << "Reading data from System " << name << ", of type " << type << " ... \n";
				ou << name << ", type " << type << " ... \n";
			}
			
			if(!type.compare("standard")) readStandard(in,dim);
			break;
		}
	}
	checkCreatedSystem();
	if (parentProc->num==rankk) ou.close();
}

// ***

MecSystem::MecSystem(string nam, map<string,MecSystem*,less<string> > mch, Process *pProc, int di)
:	name(nam),
	type("standard"),
	dim(di),
	mBodies(),
	mPoints(),
	childSysNamesG(),
	mChildSys(mch),
	parentSysName(""),
	parentProc(pProc),
	isWellDefined(false),
	pointSets(),
	mPointIniPos(),
	mPointIniVel()
{}

// ***

void MecSystem::readStandard(ifstream &in,int dim)
{
	char keyword[20];
	string auxs;
	string auxs1;
	vector<double> auxv(dim);
	vector<double> value;
	double val;
	int i;
	Joint *kkJ;
	vector<string> auxsv;
	
	while (in >> keyword && strcmp(keyword,"ENDSYSTEM") )
	{
		if(!strcmp(keyword,"POINT"))
		{
			if(dim==1)
			{
				in >> auxs >> auxv[0];
				addSystemName(auxs,name);
			}
			else if(dim==2)
			{
				in >> auxs >> auxv[0] >> auxv[1];
				addSystemName(auxs,name);
			} 
			else if(dim==3)
			{
				in >> auxs >> auxv[0] >> auxv[1] >> auxv[2];
				addSystemName(auxs,name);
			}
			pointNames[auxs] = auxv;
		}
		else if(!strcmp(keyword,"BEGINBODY"))
		{
			in >> auxs;			/** Here I read keyword NAME */
			in >> auxs;			/** Here I read the actual name */
			addSystemName(auxs,name);
			in >> auxs1;			/** Here I read keyword TYPE*/
			in >> auxs1;			/** Here I read the actual type */
			if(!auxs1.compare("flexchain")) 
			{
				mBodies[auxs] = new FlexChain(in,auxs,this);
			}
			else if(!auxs1.compare("lattice")) 
			{
				mBodies[auxs] = new Lattice(in,auxs,this);
			}
			else if(!auxs1.compare("flexfem")) 
			{
				mBodies[auxs] = new FlexFem(in,auxs,this);
			}
			else if(!auxs1.compare("particle")) 
			{
				mBodies[auxs] = new ParticleBody(in,auxs,this);				
			}
			else if(!auxs1.compare("rigidbody3d4p")) 
			{
				mBodies[auxs] = new RigidBody3d4p(in,auxs,this);
				auxsv.clear();
				auxsv.push_back(auxs+"-0");
				auxsv.push_back(auxs+"-1");
				auxsv.push_back(auxs+"-2");
				auxsv.push_back(auxs+"-3");
				mJoints[auxs+"-j"] = new RigidTetra(in,auxs+"-j",this, auxsv);
			}
			else if(!auxs1.compare("ground")) 
			{
				mBodies[auxs] = new Ground(in,auxs,this);				
			}
			else if(!auxs1.compare("thermodamper")) 
			{
				mBodies[auxs] = new ThermoDamper(in,auxs,this);
			}
			else if(!auxs1.compare("torsionfreebeam")) 
			{
				readTorsionFreeBeam(in,auxs);
			}			
			else if(!auxs1.compare("thermostickslip")) 
			{
				mBodies[auxs] = new ThermoStickSlip(in,auxs,this);
			}
            else if(!auxs1.compare("rigidbar"))
			{
                readRigidBar(in,auxs);
			}
			// .... Other types of bodies
//			else if 
			else if(parentProc->num==0) cout << "	*ERROR: Body " << auxs << " has wrong type, not added to system\n";
			
		}
		else if(!strcmp(keyword,"BEGINJOINT"))
		{
			in >> auxs;			/** Here I read keyword NAME */
			in >> auxs;			/** Here I read the actual name */
			addSystemName(auxs,name);
			in >> auxs1;			/** Here I read keyword TYPE*/
			in >> auxs1;			/** Here I read the actual type */
			if(!auxs1.compare("encastre"))
			{
				mJoints[auxs] = new Encastre(in,auxs,this);
			}
			if(!auxs1.compare("encastrepenalty"))
			{
				mJoints[auxs] = new EncastrePenalty(in,auxs,this);
			}			
			else if(!auxs1.compare("glue"))
			{
//				mJoints[auxs] = new Glue(in,auxs,this);
				kkJ = new GGlue(in,auxs,this);
			}
			/*
			else if(!auxs1.compare("vl")) 
			{
				mJoints[auxs] = new Vl(in,auxs,this);
			}
			*/
			/*
			else if(!auxs1.compare("thermodamper")) 
			{
				mJoints[auxs] = new ThermoDamper(in,auxs,this);
//				val = 300.0;
//				setTemperature(val);
			}
			*/
			else if(!auxs1.compare("rigidlink")) 
			{
				mJoints[auxs] = new RigidLink(in,auxs,this);
			}
			else if(!auxs1.compare("contact")) 
			{
				mJoints[auxs] = new Contact(in,auxs,this);
			}									
		}
		else if(!strcmp(keyword,"BEGINPOINTSET"))
		{
			in >> auxs;			/** Here I read keyword NAME */
			in >> auxs;			/** Here I read the actual name */
			addSystemName(auxs,name);
			in >> auxs1;			/** Here I read keyword X,Y,Z,XY,XZ,YZ,XYZ */
//			cout << "ausx1,size es " << auxs1.size() << "\n";
			value.clear();
			for(i=0; i<auxs1.size();++i)
			{
				in >> val;			/** Here I read the value */
				value.push_back(val);
			}
			setPointSet(auxs, auxs1, value);
			
			if(rankk==parentProc->num) cout << "Reading pointSet " << auxs << endl;
			//map<string,vector<string>,less<string> >::iterator kk;
			//int kkj;
			//for(kk=pointSets.begin(); kk!=pointSets.end(); ++kk)
			//{
//				for(kkj=0; kkj<(kk->second).size(); ++kkj)	cout << "Set " << kk->first << " contains points " << kk->second[kkj] << "\n";
			//}
			
		}		
		else if(!strcmp(keyword,"SYSTEM"))
		{
			in >> auxs;
			childSysNamesG.push_back(auxs);		/** Not checking at this stage */
		}
	}
}

// ***

void MecSystem::checkCreatedSystem()
{
	if(mBodies.size()==0 && childSysNamesG.size()==0) {
		isWellDefined = false;
		if(parentProc->num==0)  cout << "*ERROR: No bodies or other systems found in system " << name << "\n";
	} else {
		isWellDefined = true;
		if(parentProc->num==rankk)
		{
			cout << "System " << name << ", of type " << type <<", have bodies and/or systems, it looks formally well defined \n";
			cout << "--------\n";
//			cout << "----" << parentProc->num << " ----\n";
		}
	}
}

// ***

MecSystem* MecSystem::getPTSystem(string sysname)
{
	if(!name.compare(sysname) ) return this;
	else return( parentProc->getPTSystem(sysname) );
}

// ***

MecSystem* MecSystem::getPTSystemLocal(string sysname)
{
	return( (parentProc->mMecSystems[sysname]) );
}

// ***

Point* MecSystem::getPTPoint(string pointName)
{
	if(checkIfPoint(pointName)==true)
	{
		Point* aux;
		string namesys = getStringHead(pointName);
		aux = &(getPTSystem(namesys)->mPoints[pointName]);
//JC			aux = getPTSystem(namesys)->mPoints[pointName];
		return(aux);
	}
	else
	{
		cout << "*ERROR: point " <<  pointName << " does not exist \n";
		return(0); 
	}
}

// ***

Point* MecSystem::getPTPointLocal(string pointName)
{
	Point* aux;
		
	// First, check if the point was originally local but it moved (for instance, because a glue joint)
	if(mChangedPoints.count(pointName)>0)	
	{
		pointName = mChangedPoints[pointName];
	}

	string namesys = getStringHead(pointName);
	
	/** If the point is from other process, look at the root system */
	if(checkIfSystem(namesys)==false)
	{
		aux = &((parentProc->rootSys)->mPoints[pointName]);
//JC			aux = parentProc->rootSys->mPoints[pointName];
	}
	else
	{
		aux = &(getPTSystemLocal(namesys)->mPoints[pointName]);
//JC			aux = getPTSystemLocal(namesys)->mPoints[pointName];
	}
	return(aux);
}

// ***

void MecSystem::setParenthood()
{
	int i;
	string auxs;
	
/** Identify the child systems that belong to this process and put the parent name if appropriate */
	for(i=0; i<childSysNamesG.size(); i++)
	{
		auxs = childSysNamesG[i];
		if( (parentProc->mMecSystems).count(auxs) )
		{
			mChildSys[auxs] = parentProc->mMecSystems[auxs];
			(mChildSys[auxs])->parentSysName = name;
		}
	}
}

// ***

void MecSystem::createPoint(const Point* ppo)
//JC void MecSystem::createPoint(Point* ppo)
{
	mPoints[ppo->name] = *ppo;
//JC	mPoints[ppo->name] = ppo;
}

// ***

void MecSystem::createPointAtRoot(const Point* ppo)
//JC void MecSystem::createPointAtRoot(Point* ppo)
{
	parentProc->rootSys->createPoint(ppo);
}

// ***

void MecSystem::addPointName(const string & name, const vector<double>& coords)
{
	pointNames[name] = coords;
}

// ***

void MecSystem::createPoints()
{
    map<string,MecSystem*,less<string> >::iterator j;
	map<string,string,less<string> >::iterator i;
	string auxs;
//	string oldName;
//	string newName;
/*	
	// Create points in the root system, which are copies of points from other processes
	if(!name.compare("root"))
	{
		cout << "Estoy en root, y mChangedPoints.size es  " << mChangedPoints.size() << endl;
		for(i=mChangedPoints.begin(); i!=mChangedPoints.end(); i++)
		{
//			oldName = i->first;
//			newName = mChangedPoints[i->first];
			mPoints[ i->second ] = *(new Point(i->second, pointNames[i->second]));
			//JC			mPoints[auxs] = new Point(auxs,k->second);
			cout << "Creado punto " << i->second << endl;
		}		
	}
*/	
    /** Recursive creation */
    for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->createPoints();
	
	map<string, vector<double>, less<string> >::iterator k;
		
//	if (parentProc->num == rankk)	cout << "Creando puntos en system " << name << endl;
	for(k=pointNames.begin(); k!=pointNames.end(); k++)
	{
//		if (parentProc->num == rankk)	cout << "creando punto " << k->first << " en system " << name << endl;
		auxs = k->first;
		mPoints[auxs] = *(new Point(auxs,k->second));
//JC			mPoints[auxs] = new Point(auxs,k->second);
	}
	// OJO!!!  Esto puede dar problemas.  Cambiar mejor a map de punteros a puntos
}


// ***

void MecSystem::setPPointersCompBodies()
{
    map<string,Body*,less<string> >::iterator i;
    map<string,MecSystem*,less<string> >::iterator j;
	
	/** Recursive setting */
    for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->setPPointersCompBodies();

	/** Request the bodies to set the pointers to points in every component */
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		i->second->setPTPointsComp();
	}
}

// ***

void MecSystem::setPPointersLoadBodies()
{
    map<string,Body*,less<string> >::iterator i;
    map<string,MecSystem*,less<string> >::iterator j;
	
	/** Recursive setting */
    for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->setPPointersLoadBodies();

	/** Request the bodies to set the pointers to points in every component */
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		(*((*i).second)).setPTPointsLoad();
	}
}

// ***
/*
void MecSystem::setPPointersCompJoints()
{
    map<string,Joint*,less<string> >::iterator i;
    map<string,MecSystem*,less<string> >::iterator j;
	
	// Recursive setting
    for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->setPPointersCompJoints();
	
	// Request the joints to set the pointers to points in every component
	for(i=mJoints.begin(); i!=mJoints.end(); i++)
	{
		(*((*i).second)).setPTPointsComp();
	}
}

// ***

void MecSystem::setPPointersLoadJoints()
{
    map<string,Joint*,less<string> >::iterator i;
    map<string,MecSystem*,less<string> >::iterator j;
	
	// Recursive setting
    for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->setPPointersLoadJoints();
	
	// Request the joints to set the pointers to points in every component
	for(i=mJoints.begin(); i!=mJoints.end(); i++)
	{
		(*((*i).second)).setPTPointsLoad();
	}
}
*/
// ***

void MecSystem::setPPointersConsJoints()
{
    map<string,Joint*,less<string> >::iterator i;
    map<string,MecSystem*,less<string> >::iterator j;
	
	/** Recursive setting */
    for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->setPPointersConsJoints();
	
	/** Request the joints to set the pointers to points in every component */
	for(i=mJoints.begin(); i!=mJoints.end(); i++)
	{
		(*((*i).second)).setPTPointsCons();
	}
}

// ***

void MecSystem::setPPointSurface()
{
    map<string,Body*,less<string> >::iterator i;
    map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive setting */
    for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->setPPointSurface();
	
	/** Request the bodies to set the pointers to points in every component */
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		i->second->setPPointSurface();
	}
	/** Request the joints to set the pointers to points in every component */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		k->second->setPPointSurface();
	}	
}

// ***

bool MecSystem::checkIfPointProcess(string st)
{
	/** First check by the system name */
	if( (parentProc->checkIfSystem(getStringHead(st)))==true ) return(true);
	else	/** next, check if it is at the root system */
	if( checkIfPointRootSystem(st)==true) return(true);
	return(false);
}

// ***

bool MecSystem::checkIfPointSystem(string st)
{	
	if(pointNames.count(st)) return(true); else return(false);
}

// ***

bool MecSystem::checkIfPointRootSystem(string st)
{
	return( parentProc->rootSys->checkIfPointSystem(st) );
}

// ***

bool MecSystem::checkIfSystem(string st)
{
	return(parentProc->checkIfSystem(st));
}

// ***

bool MecSystem::checkIfPoint(string st)
{
	string sysName = getStringHead(st);
	MecSystem* psys=parentProc->getPTSystem(sysName);
	return(psys->checkIfPointSystem(st));
}

// ***

bool MecSystem::checkIfPointChanged(string st)
{
	if(mChangedPoints.count(st)>0) return(true); else return(false);
}

// ***

void MecSystem::addRemoteSharedPoint(string pointName, int numProc)
{
	parentProc->addRemoteSharedPoint(pointName,numProc);
}

// ***

int MecSystem::getProcNum()
{
	return(parentProc->num);
}

// ***

int MecSystem::getRankNum()
{
	return(parentProc->lrank);
}

// ***

void MecSystem::setDofsBodies()
{
	map<string,Body*,less<string> >::iterator i;
	map<string,MecSystem*,less<string> >::iterator k;

    for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->setDofsBodies();	/** Recursion */
	
	if(mBodies.size()!=0 )		
	{
		for(i=mBodies.begin(); i!=mBodies.end(); i++)	((*i).second)->setDofs();
	}
}

// ***

void MecSystem::setDofsJoints()
{
	map<string,Joint*,less<string> >::iterator j;
	map<string,MecSystem*,less<string> >::iterator k;
	
    for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->setDofsJoints();	/** Recursion */

	if(mJoints.size()!=0 )		
	{
		for(j=mJoints.begin(); j!=mJoints.end(); j++)	((*j).second)->setDofs();
	}	
}

// ***

int MecSystem::countDofs()
{
/*
	//Recursive	
	int numDofs=0;
	map<string,Point,less<string> >::iterator j;
//JC	map<string,Point*,less<string> >::iterator j;
	map<string,MecSystem*,less<string> >::iterator k;

    for(k=mChildSys.begin(); k!=mChildSys.end(); k++) numDofs+= ( ((*k).second)->countDofs() );
	
	// Loop on its points
	for(j=mPoints.begin(); j!=mPoints.end(); j++)	numDofs+= ((*j).second).countDofs();
//JC		for(j=mPoints.begin(); j!=mPoints.end(); j++)	numDofs+= j->second->countDofs();
	return(numDofs);
*/
}

// ***

int MecSystem::setTransNumbers(int num)
{
	//Recursive	
	map<string,Point,less<string> >::iterator j;
//JC	map<string,Point*,less<string> >::iterator j;
	map<string,MecSystem*,less<string> >::iterator k;

    for(k=mChildSys.begin(); k!=mChildSys.end(); k++) num = ( ((*k).second)->setTransNumbers(num) );	/** Recursion */
	// Loop on its points
	for(j=mPoints.begin(); j!=mPoints.end(); j++)	num = ((*j).second).setTransNumbers(num);
//JC		for(j=mPoints.begin(); j!=mPoints.end(); j++)	num = j->second->setTransNumbers(num);
	return(num);
	
	// Account for component dofs. TO BE DONE
}

// ***

int MecSystem::setRotNumbers(int num)
{
	//Recursive	
	map<string,Point,less<string> >::iterator j;
//JC	map<string,Point*,less<string> >::iterator j;
	map<string,MecSystem*,less<string> >::iterator k;

    for(k=mChildSys.begin(); k!=mChildSys.end(); k++) num = ( ((*k).second)->setRotNumbers(num) );	/** Recursion */
	// Loop on its points
	for(j=mPoints.begin(); j!=mPoints.end(); j++)	num = ((*j).second).setRotNumbers(num);
//JC		for(j=mPoints.begin(); j!=mPoints.end(); j++)	num = j->second->setRotNumbers(num);
	return(num);
	
	// Account for component dofs. TO BE DONE
}

// ***

int MecSystem::setTRNumbers(int num)
{
	//Recursive	
	map<string,Point,less<string> >::iterator j;
	map<string,MecSystem*,less<string> >::iterator k;

    for(k=mChildSys.begin(); k!=mChildSys.end(); k++) num = ( ((*k).second)->setTRNumbers(num) );	/** Recursion */
	
	// Loop on its points
	for(j=mPoints.begin(); j!=mPoints.end(); j++)	num = ((*j).second).setTRNumbers(num);
	
	// Loop on bodies, looking for component dofs
	map<string,Body*,less<string> >::iterator itb;
	map<string,Joint*,less<string> >::iterator itj;

//cout << "num A = " << num << endl;

	for(itb=mBodies.begin(); itb!=mBodies.end(); itb++)
	{
		num = itb->second->setTRNumbers(num);
//		cout << "UUUUU " << num << endl;
	}
/*
	cout << "num B = " << num << endl;
	for(itj=mJoints.begin(); itj!=mJoints.end(); itj++)
	{
		itj->second->setTRNumbers(num);
//		cout << "UUUUU " << num << endl;
	}
*/
//cout << "num C = " << num << endl;
//cout << "mJoints = " << mJoints.size() << endl;
	
	return(num);
}

// ***

void MecSystem::setLocalGlob()
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator j;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->setLocalGlob();	/** Recursion */

	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		((*i).second)->setLocalGlob();
	}
	for(j=mJoints.begin(); j!=mJoints.end(); j++)
	{
		((*j).second)->setLocalGlob();
	}
}

// ***

void MecSystem::initializeGlobalDofsBodies(vector<double> &q, vector<double> &dq, double & time)
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->initializeGlobalDofsBodies(q, dq, time);	/** Recursion */

	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		((*i).second)->initializeGlobalDofs(q, dq, time);
	}
}

// ***

void MecSystem::initializeGlobalDofsJoints(vector<double> &q, vector<double> &dq, double & time)
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Joint*,less<string> >::iterator j;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->initializeGlobalDofsJoints(q, dq, time);	/** Recursion */
	
	for(j=mJoints.begin(); j!=mJoints.end(); j++)
	{
		((*j).second)->initializeGlobalDofs(q, dq, time);
	}	
}

// ***

void MecSystem::setMatAtComp()
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator j;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->setMatAtComp();	/** Recursion */

	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		((*i).second)->setMatAtComp();
	}

	// There are joints that need paramaters (rigidchain for torsionFreeBemas for instance)
	for(j=mJoints.begin(); j!=mJoints.end(); j++)
	{
		((*j).second)->setParamAtCons();
	}
}

// ***

Material* MecSystem::getPTMaterial(const string& matName)
{
	return parentProc->mMaterials[matName];
}

// ***

void MecSystem::fillResidual(vector<double>& res_in, vector<double>& q_in, vector<double>& dq_in, const double& time, const int& modeFlag, const double& grav, const int& gravDir, const double& DeltaT)
{

	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator j;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->fillResidual(res_in, q_in, dq_in,time, modeFlag, grav, gravDir, DeltaT);	/** Recursion */

	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		((*i).second)->fillResidual(res_in, q_in, dq_in, time, modeFlag, grav, gravDir, DeltaT);
	}
	for(j=mJoints.begin(); j!=mJoints.end(); j++)
	{
		((*j).second)->fillResidual(res_in, q_in, dq_in, time, modeFlag, grav, gravDir,DeltaT);
	}	
}

// ***

void MecSystem::fillResidual(	vector<double>& residue,
								const vector<double>& q,
								const vector<double>& qdot,
								const vector<double>& qddot,
								const double& time,
								const int& modeFlag,
								const double& grav, const int& gravDir,
//								const int& integratorType, const vector<double>& qConv, const vector<double>& dqConv,
 								const int& integratorType,
								const double& DeltaT
							)
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator j;
	
//JC	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->fillResidual(residue, q, qdot, qddot, time, modeFlag, grav, gravDir, integratorType, qConv, dqConv, DeltaT);
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->fillResidual(residue, q, qdot, qddot, time, modeFlag, grav, gravDir, integratorType, DeltaT);

	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
//JC		((*i).second)->fillResidual(residue, q, qdot, qddot, time, modeFlag, grav, gravDir, integratorType, qConv, dqConv, DeltaT);
		((*i).second)->fillResidual(residue, q, qdot, qddot, time, modeFlag, grav, gravDir, integratorType, DeltaT);
	}
	for(j=mJoints.begin(); j!=mJoints.end(); j++)
	{
//JC		((*j).second)->fillResidual(residue, q, qdot, qddot, time, modeFlag, grav, gravDir, integratorType, qConv, dqConv, DeltaT);
		((*j).second)->fillResidual(residue, q, qdot, qddot, time, modeFlag, grav, gravDir, integratorType, DeltaT);
	}	
}

// ***

void MecSystem::fillResFi(	vector<double>& residue,
							 const vector<double>& q,
							 const vector<double>& qddot,
							 const int& modeFlag
							 )
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator j;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->fillResFi(residue, q, qddot, modeFlag);
	
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		((*i).second)->fillResFi(residue, q, qddot, modeFlag);
	}
	for(j=mJoints.begin(); j!=mJoints.end(); j++)
	{
		((*j).second)->fillResFi(residue, q, qddot, modeFlag);
	}	
}

// ***

void MecSystem::fillResFq(	vector<double>& residue,
							 const vector<double>& q,
							 const vector<double>& qdot,
							 const int& modeFlag,
//							 const int& integratorType, const vector<double>& qConv, const vector<double>& dqConv,
  							 const int& integratorType,
							 const double& DeltaT
							 )
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator j;
	
//JC	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->fillResFq(residue, q, qdot, modeFlag, integratorType, qConv, dqConv, DeltaT);
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->fillResFq(residue, q, qdot, modeFlag, integratorType, DeltaT);
	
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
//JC		((*i).second)->fillResFq(residue, q, qdot, modeFlag, integratorType, qConv, dqConv, DeltaT);
		((*i).second)->fillResFq(residue, q, qdot, modeFlag, integratorType, DeltaT);
	}
	for(j=mJoints.begin(); j!=mJoints.end(); j++)
	{
//JC		((*j).second)->fillResFq(residue, q, qdot, modeFlag, integratorType, qConv, dqConv, DeltaT);
		((*j).second)->fillResFq(residue, q, qdot, modeFlag, integratorType, DeltaT);
	}		
}

// ***

void MecSystem::fillResFt(	vector<double>& residue,
							const double& time,
							const int& modeFlag,
							const double& grav, const int& gravDir
							)
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator j;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->fillResFt(residue, time, modeFlag, grav, gravDir);
	
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		((*i).second)->fillResFt(residue, time, modeFlag, grav, gravDir);
	}
	for(j=mJoints.begin(); j!=mJoints.end(); j++)
	{
		((*j).second)->fillResFt(residue, time, modeFlag, grav, gravDir);
	}			
}
	
// ***
/*
void MecSystem::fillTangent(sparseMatrix& jac_in, vector<double>& q_in, const double& time, const int& modeFlag)
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator j;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->fillTangent(jac_in, q_in, time, modeFlag);

	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		((*i).second)->fillTangent(jac_in, q_in, time, modeFlag);
	}
	for(j=mJoints.begin(); j!=mJoints.end(); j++)
	{
		((*j).second)->fillTangent(jac_in, q_in, time, modeFlag);
	}		
}
*/
// ***

void MecSystem::fillTangent(sparseMatrix& Kint, vector<vector<double> >& Kbb, vector<vector<double> >& Kib, vector<vector<double> >& Kbi,
							vector<int>& globalLocalI, vector<int>& globalLocalB,
							vector<double>& q_in, vector<double>& dq_in,
							const double& time, const int& modeFlag)
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator j;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->fillTangent(Kint, Kbb, Kib, Kbi, globalLocalI, globalLocalB, q_in, dq_in, time, modeFlag);
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		((*i).second)->fillTangent(Kint, Kbb, Kib, Kbi, globalLocalI, globalLocalB, q_in, dq_in, time, modeFlag);
	}
	for(j=mJoints.begin(); j!=mJoints.end(); j++)
	{
		((*j).second)->fillTangent(Kint, Kbb, Kib, Kbi, globalLocalI, globalLocalB, q_in, dq_in, time, modeFlag);
	}			
}

// ***
/*
void MecSystem::fillTangent	(	sparseMatrix& tangent,
								const vector<double>& q,
								const vector<double>& qdot,
								const vector<double>& qddot,
								double partial_qdot,
								double partial_qddot,
								const double& time,
								const int& modeFlag,
 								const int& integratorType,
								const double& DeltaT
							)
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator j;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->fillTangent(tangent, q, qdot, qddot, partial_qdot, partial_qddot, time, modeFlag, integratorType, DeltaT);

	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		((*i).second)->fillTangent(tangent, q, qdot, qddot, partial_qdot, partial_qddot, time, modeFlag, integratorType, DeltaT);
	}
	for(j=mJoints.begin(); j!=mJoints.end(); j++)
	{
		((*j).second)->fillTangent(tangent, q, qdot, qddot, partial_qdot, partial_qddot, time, modeFlag, integratorType, DeltaT);
	}		
}
*/

// ***

void MecSystem::fillTangent	(sparseMatrix& Kint, vector<vector<double> >& Kbb, vector<vector<double> >& Kib, vector<vector<double> >& Kbi,
							 vector<int>& globalLocalI, vector<int>& globalLocalB,	
							 const vector<double>& q,
							 const vector<double>& qdot,
							 const vector<double>& qddot,
							 double partial_qdot,
							 double partial_qddot,
							 const double& time,
							 const int& modeFlag,
							 //								const int& integratorType, const vector<double>& qConv, const vector<double>& dqConv,
							 const int& integratorType,
							 const double& DeltaT
							 )
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator j;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->fillTangent(Kint, Kbb, Kib, Kbi, globalLocalI, globalLocalB, q, qdot, qddot, partial_qdot, partial_qddot, time, modeFlag, integratorType, DeltaT);	/** Recursion */
	
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		((*i).second)->fillTangent(Kint, Kbb, Kib, Kbi, globalLocalI, globalLocalB, q, qdot, qddot, partial_qdot, partial_qddot, time, modeFlag, integratorType, DeltaT);
	}
	for(j=mJoints.begin(); j!=mJoints.end(); j++)
	{
		((*j).second)->fillTangent(Kint, Kbb, Kib, Kbi, globalLocalI, globalLocalB, q, qdot, qddot, partial_qdot, partial_qddot, time, modeFlag, integratorType, DeltaT);
	}		
}

// ***

void MecSystem::fillMassMatrix(sparseMatrix& Ma, const vector<double>& q, double time)
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator j;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->fillMassMatrix(Ma,q,time);	/** Recursion */

	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		((*i).second)->fillMassMatrix(Ma,q,time);
	}
	for(j=mJoints.begin(); j!=mJoints.end(); j++)
	{
		((*j).second)->fillMassMatrix(Ma,q,time);
	}				
}

// ***

void MecSystem::fillMassMatrix(vector<double>& vMa, const vector<double>& q, double time)
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator j;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->fillMassMatrix(vMa,q,time);	/** Recursion */
	
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		((*i).second)->fillMassMatrix(vMa,q,time);
	}
	for(j=mJoints.begin(); j!=mJoints.end(); j++)
	{
		((*j).second)->fillMassMatrix(vMa,q,time);
	}				
}

// ***
void MecSystem::fillForce(vector<double>& f,
						const vector<double>& q,
						const vector<double>& qdot,
						const double& time,
						const double& grav, const int& gravDir
					)
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator j;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->fillForce(f, q, qdot, time, grav, gravDir);	/** Recursion */

	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		((*i).second)->fillForce(f, q, qdot, time, grav, gravDir);
	}
	for(j=mJoints.begin(); j!=mJoints.end(); j++)
	{
		((*j).second)->fillForce(f, q, qdot, time, grav, gravDir);
	}					
}

// ***

void MecSystem::reportContent()
{
	map<string,Body*,less<string> >::iterator i;
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Point,less<string> >::iterator k;
	map<string,Joint*,less<string> >::iterator l;
	
	if(mBodies.size()!=0 || mJoints.size()!=0 ||mChildSys.size()!=0 || mPoints.size()!=0 )
	{
		cout <<	"	System " << name << " (" << type << ") contains \n";
		/** Reports points */
		for(k=mPoints.begin(); k!=mPoints.end(); k++)
		{
			cout << "		Point " << (k->second).name << " = [";
			for(int n=0;n<dim-1;n++) cout << (k->second).r[n] << ",";
			cout << (k->second).r[dim-1] << "] \n";
		}
		/** Reports bodies */
		for(i=mBodies.begin(); i!=mBodies.end(); i++)
		{
			cout << "		Body " << ((*i).second)->name << " contains \n";
			/** Reports components within the body */
			((*i).second)->reportContent();
		}
		/** Reports joints */
		for(l=mJoints.begin(); l!=mJoints.end(); l++)
		{
			cout << "		Joint " << ((*l).second)->name << " contains \n";
			/** Reports components within the joint */
			((*l).second)->reportContent();
		}						
		/** Reports systems */
		for(j=mChildSys.begin(); j!=mChildSys.end(); j++)
		{
			cout << "		System " << (j->second)->name << "\n";
		}
	}
	else
	{
		cout << "*WARNING: System " << name << " does not have any valid body or child system in this process \n";
	}

	/** Recursive report */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->reportContent();
}

// ***

void MecSystem::reportPointsDofs()
{
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Point,less<string> >::iterator k;
	int i;
	int nt;
	int nr;
	
	for(k=mPoints.begin(); k!=mPoints.end(); k++)
	{
		cout << (k->second).name << "\t \t ";
		cout << parentProc->num << "\t \t ";
		nt = (k->second).globalTNum.size();
		nr = (k->second).globalRNum.size();
		for(i=0;i<nt;i++) cout << (k->second).globalTNum[i] << "\t";
		for(i=0;i<nr;i++) cout << (k->second).globalRNum[i] << "\t";
		cout << "\n";
	}
		/** Recursive report */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->reportPointsDofs();
}

// ***

void MecSystem::reportPointsDofsFile()
{
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Point,less<string> >::iterator k;
	int i;
	int nt;
	int nr;
	
	/** Recursive report */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->reportPointsDofsFile();

//	if (!name.compare("e1") && parentProc->num==rankk)
//	{
//		cout << "EN sistema " << name << " de proceso " << rankk << " mPoints size es " << mPoints.size() <<  endl;
//	}
//	int kk=0;
	for(k=mPoints.begin(); k!=mPoints.end(); k++)
	{
//		if (!name.compare("e1") && parentProc->num==rankk)
//		{
//			cout << "kk " << kk << " EN sistema " << name << " de proceso " << rankk << " mPoint  es " << k->first <<  endl;
//		}
//		kk+=1;
		
		parentProc->writeToFile0( (k->second).name );
		parentProc->writeToFile0( "\t " );
		nt = (k->second).globalTNum.size();
		nr = (k->second).globalRNum.size();
		for(i=0;i<nt;i++) 
		{
			parentProc->writeToFile( (k->second).globalTNum[i],"\t");
		}
		for(i=0;i<nr;i++)
		{
			parentProc->writeToFile( (k->second).globalRNum[i],"\t");
		}
		parentProc->writeToFile0( "\n" );
	}
}

// ***

void MecSystem::writeFileGraph(ofstream& fGraph, string& gFormat, string& gPath, string& aux, bool& plotDispFlag, bool& plotProcFlag, bool& plotStressFlag)
{
	map<string,Body*,less<string> >::iterator i;
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->writeFileGraph(fGraph,gFormat,gPath,aux,plotDispFlag,plotProcFlag,plotStressFlag);

	/** Request the bodies to provide graphical output */
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		(*((*i).second)).writeFileGraph(fGraph,gFormat,gPath,aux,plotDispFlag,plotProcFlag,plotStressFlag);
	}
	/** Request the joints to provide graphical output */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		(*((*k).second)).writeFileGraph(fGraph,gFormat,gPath,aux);
	}	
}

// ***

double MecSystem::getDof(int num)
{
	return parentProc->q[num];
}

// ***

void MecSystem::setDofsRoot()
{
	
	map<string,Point,less<string> >::iterator i;
	Point* remotePoint;
	
	for (i=mPoints.begin(); i!=mPoints.end(); ++i)
	{
		if ( (i->second).countDofs()==0 )
		{
			remotePoint = getPTPoint(i->first);
			(i->second).setDofT( remotePoint->getDofT() );
			(i->second).setDofR( remotePoint->getDofR() );
		}
	}
}

// ***

vector<string> MecSystem::getPointNames()
{
	vector<string> v;
	map<string,Point,less<string> >::iterator i;
	for (i=mPoints.begin(); i!=mPoints.end(); ++i)
	{
		v.push_back(i->first);
	}
	return(v);
}

// ***

void MecSystem::updateDofs(vector<double> & q, vector<double> & deltaq)
{
	/** Recursive */
	map<string,MecSystem*,less<string> >::iterator j;
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->updateDofs(q,deltaq);

	map<string,Point,less<string> >::iterator i;
	for (i=mPoints.begin(); i!=mPoints.end(); ++i)
	{
		(i->second).updateDofs(q,deltaq);
	}
	
	// Update component dofs
	
	map<string,Body*,less<string> >::iterator iteb;
	map<string,Joint*,less<string> >::iterator itej;

	/** Request the bodies */
	for(iteb=mBodies.begin(); iteb!=mBodies.end(); iteb++)
	{
		iteb->second->updateDofs(q,deltaq);
	}

	/** Request the joints */
	for(itej=mJoints.begin(); itej!=mJoints.end(); itej++)
	{
		itej->second->updateDofs(q,deltaq);
	}	
}

// ***

void MecSystem::setSparsePattern(vector<vector<int> >& tempIa)
{
	map<string,Body*,less<string> >::iterator i;
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->setSparsePattern(tempIa);
	
	/** Request the bodies */
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		(*((*i).second)).setSparsePattern(tempIa);
	}

	/** Request the joints */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		(*((*k).second)).setSparsePattern(tempIa);
	}	
}

// ***

void MecSystem::setSparsePatternInertia(vector<vector<int> >& tempIa)
{
	map<string,Body*,less<string> >::iterator i;
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->setSparsePatternInertia(tempIa);
	
	/** Request the bodies */
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		(*((*i).second)).setSparsePatternInertia(tempIa);
	}
	
	/** Request the joints */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		(*((*k).second)).setSparsePatternInertia(tempIa);
	}	
}

// ***

void MecSystem::changePointName(string oldName, string newName)
{
	//cout << "En System " << name << endl;
	//cout << "oldname is " << oldName << endl;
	//cout << "newname is " << newName << endl;
	
	// Only name is changed.  Coordinates are unchanged
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator k;
	vector<double> coords;
	
	// First, I check that the point name has not already changed (by a glue joint for instance)
	if (checkIfPointSystem(oldName)==0)	// name does not exist, it has been changed
	{
		oldName = mChangedPoints[oldName];
	}
	
//	if (rankk==parentProc->num)  cout << "Changing point name. Old is " << oldName << "; new is " << newName << "\n";
	coords = pointNames[oldName];
	
//	cout << "Coordenadas antiguo " << pointNames[oldName][0] << "\t" << pointNames[oldName][1] << "\t" << pointNames[oldName][2] << "\n";		
	pointNames.erase(oldName);
	
	pointNames[newName] = coords;	
	
//	cout << "Coordenadas nuevo " << pointNames[newName][0] << "\t" << pointNames[newName][1] << "\t" << pointNames[newName][2] << "\n";		
	mChangedPoints[oldName]=newName;
}

// ***

char* MecSystem::getOutGraphFile()
{
	return parentProc->outGraphFile;
}

// ***

void MecSystem::initializeCompRes(const bool& st, const double& grav, const int& gravDir, const bool augLagFlag)
{
	map<string,Body*,less<string> >::iterator i;
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->initializeCompRes(st,grav,gravDir,augLagFlag);
	
	/** Request the bodies */
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		(*((*i).second)).initializeCompRes(st,grav,gravDir);
	}
	/** Request the joints */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		(*((*k).second)).initializeCompRes(st,grav,gravDir,augLagFlag);
	}		
}

// ***

//void MecSystem::initializeCompTang(const bool& st, const bool& lmf)
void MecSystem::initializeCompTang(const bool& st)
{
	
	map<string,Body*,less<string> >::iterator i;
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
//	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->initializeCompTang(st,lmf);
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->initializeCompTang(st);
	
	/** Request the bodies */
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
//		(*((*i).second)).initializeCompTang(st,lmf);
		(*((*i).second)).initializeCompTang(st);
	}
	/** Request the joints */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
//		(*((*k).second)).initializeCompTang(st,lmf);
		(*((*k).second)).initializeCompTang(st);
	}		
}

// ***

void MecSystem::setPointSet(string& nameSet, string& coor, vector<double>& value)
{
	// I have to search in pointNames
	map<string, vector<double>, less<string> >::iterator i;
	int j, k;
	vector<int> icoor;
	
	if(!coor.compare("1"))	icoor.push_back(0);
	else if(!coor.compare("2"))	icoor.push_back(1);
	else if(!coor.compare("3"))	icoor.push_back(2);
	else if(!coor.compare("12"))
	{
		icoor.push_back(0);
		icoor.push_back(1);
	}
	else if(!coor.compare("13"))
	{
		icoor.push_back(0);
		icoor.push_back(2);
	}
	else if(!coor.compare("23"))
	{
		icoor.push_back(1);
		icoor.push_back(2);
	}
	else if(!coor.compare("123"))
	{
		icoor.push_back(0);
		icoor.push_back(1);
		icoor.push_back(2);
	}
		
	for(i=pointNames.begin();i!=pointNames.end();++i)
	{
		if(icoor.size()==1)
		{
			if(i->second[icoor[0]]==value[0])
			{
				pointSets[nameSet].push_back(i->first);
			}
		}
		if(icoor.size()==2)
		{
			if( i->second[icoor[0]]==value[0] && i->second[icoor[1]]==value[1])
			{
				pointSets[nameSet].push_back(i->first);
			}
		}
		if(icoor.size()==3)
		{
			if(i->second[icoor[0]]==value[0] && i->second[icoor[1]]==value[1] && i->second[icoor[2]]==value[2])
			{
			   pointSets[nameSet].push_back(i->first);
			}
		}
	}
}

// ***
	   
vector<double> MecSystem::getCoordsPoint(string& st)
{
//	cout << "Accediendo a las coordenadas del punto " << st << endl;
	
		// First, check if the name exist
	   if (checkIfPointSystem(st)==1) return pointNames[st];
	   else	if(mChangedPoints.count(st)==1) // It may have changed. For instance, because a glue joint
	   {
				return pointNames[mChangedPoints[st]];
	   }
	   else
	   {
		   cout << "	*ERROR in getting point coordinates \n";
	   }
}

// ***

void MecSystem::calcEner(double& Ekin, double& Epot, double& Econ, double& Eext, double& Eamb, double& Dint, double& grav, int& gravDir, vector<double>& q, vector<double>& dq)
{
	map<string,Body*,less<string> >::iterator i;
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->calcEner(Ekin, Epot, Econ, Eext, Eamb, Dint, grav, gravDir, q, dq);
	
	/** Request the bodies */
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		(*((*i).second)).calcEner(Ekin, Epot, Econ, Eext, Eamb, Dint, grav, gravDir, q, dq);
	}
	/** Request the joints */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		(*((*k).second)).calcEner(Ekin, Epot, Econ, Eext, Eamb, Dint, grav, gravDir, q, dq);
	}			
}

// ***

void MecSystem::calcEntropy(double& Entropy, double& EntropyAmb)
{
	map<string,Body*,less<string> >::iterator i;
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->calcEntropy(Entropy, EntropyAmb);
	
	/** Request the bodies */
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		(*((*i).second)).calcEntropy(Entropy,EntropyAmb);
	}
	/** Request the joints */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		(*((*k).second)).calcEntropy(Entropy,EntropyAmb);
	}			
}

// ***

void MecSystem::calcAM(vector<double>& am, vector<double>& q_in, vector<double>& dq_in)
{
	map<string,Body*,less<string> >::iterator i;
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->calcAM(am, q_in, dq_in);
	
	/** Request the bodies */
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		(*((*i).second)).calcAM(am, q_in, dq_in);
	}
	/** Request the joints */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		(*((*k).second)).calcAM(am, q_in, dq_in);
	}			
}

// ***

void MecSystem::calcCon(vector<double> & vcon, vector<double>& q, vector<double>& dq, vector<double>& ddq, double& time)
{
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;

	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->calcCon(vcon, q, dq, ddq, time);
	
	/** Request the joints */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		(*((*k).second)).calcCon(vcon, q, dq, ddq, time);
	}				
}

// ***

void MecSystem::setIniPosPoint(string& iniPointName, vector<double>& iniValues)
{
	string kkst;
	vector<double> kkv;
	
	kkst = iniPointName;
	kkv = iniValues;

//	if (checkIfPointChanged(kkst)==true)	kkst=mChangedPoints[kkst];
	
	mPointIniPos[kkst] = kkv;
}

// ***

void MecSystem::setIniPosPointSet(string& iniPointSetName, vector<double>& iniValues)
{
	int i;
	string kkst;
	vector<double> kkv;
	
	for(i=0; i!=pointSets[iniPointSetName].size(); i++)
	{
		kkst = pointSets[iniPointSetName][i];
		kkv = iniValues;
		mPointIniPos[kkst] = kkv;
	}
}

// ***

void MecSystem::setIniVelPoint(string& iniPointName, vector<double>& iniValues)
{
//	cout << "iniPointName " << iniPointName << endl;
//	cout << "iniValues " << iniValues[0] << "\t" << iniValues[1] << endl;
	string kkst;
	vector<double> kkv;
	
	kkst = iniPointName;
	kkv = iniValues;
/*
	if (checkIfPointChanged(kkst)==true)
	{
		kkst=mChangedPoints[kkst];
		cout << "Ojo, ha cambiado a " << kkst << endl;
		if(getStringHead(kkst)!=name)	// The point is at root system
		{
			
		}
	}
*/
	mPointIniVel[kkst] = kkv;
}

// ***

void MecSystem::setIniVelPointSet(string& iniPointSetName, vector<double>& iniValues)
{
	int i;
	string kkst;
	vector<double> kkv;
	
	for(i=0; i!=pointSets[iniPointSetName].size(); i++)
	{
		kkst = pointSets[iniPointSetName][i];
		kkv = iniValues;
		mPointIniVel[kkst] = kkv;
		//cout << "kkst = " << kkst << endl;
	}
}

// ***

void MecSystem::setIniVelBody(string& bodyName, string& refPointName, vector<double>& v0, vector<double>& Omega)
{
	vector<string> vPointNames;
	vector<string>::iterator p;
	vector<double> refCoords;
	
	mBodies[bodyName]->getPointNames(vPointNames);
	// vPointNames may contain repeated names
//	cout << "vPointNames are = " << endl;
//	printSTLVector(vPointNames);

	// Order and remove duplicated names
	sort(vPointNames.begin(),vPointNames.end());
	p=unique(vPointNames.begin(),vPointNames.end());
	vPointNames.erase(p,vPointNames.end());

	// Get the coordinates of the reference point
	// The point may belog to other system and process
	refCoords = ( getPTSystem( getStringHead(refPointName) ) )->getCoordsPoint(refPointName);
	//cout << "refCoords = " << endl;
	//printSTLVector(refCoords);
	
	vector<double> vCoords(refCoords.size());
	
	// Loop on the points
	for (p=vPointNames.begin(); p!=vPointNames.end(); ++p)
	{
		// Get point coordinates
		vCoords = getCoordsPoint(*p);
		//cout << "pointName = " << *p << endl;
		//cout << "vCoords = " << endl;
		//printSTLVector(vCoords);
		substractSTLVectors(vCoords,refCoords);	// relPos stored in vCoords now
		// v = v0 + Omega x relPos
		crossProdSTLVectors(Omega,vCoords);		// (Omega x resPos) is in vCoords
		addSTLVectors(vCoords,v0);				// Velocity is in vCoords
		//cout << "velocity = " << endl;
		//printSTLVector(vCoords);
		mPointIniVel[*p] = vCoords;
	}
}

// ***

void MecSystem::setInitialConditions(vector<double>& q, vector<double>& dq)
{
	
	map<string,MecSystem*,less<string> >::iterator j;
	
    /** Recursive creation */
    for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->setInitialConditions(q,dq);
	
	map<string, vector<double>, less<string> >::iterator it;
	vector<int> gn;
	int i;
	
	// Position
	for(it=mPointIniPos.begin(); it!=mPointIniPos.end(); ++it)
	{
		gn = mPoints[it->first].getGlobalNumbers();
		for(i=0; i<gn.size(); ++i)
		{
			q[ gn[i] ] = (it->second)[i];
		}
		gn.clear();
	}
	
	// Velocity
	for(it=mPointIniVel.begin(); it!=mPointIniVel.end(); ++it)
	{
		gn = mPoints[it->first].getGlobalNumbers();
		for(i=0; i<gn.size(); ++i)
		{
			dq[ gn[i] ] = (it->second)[i];
		}
		gn.clear();
	}
	
	// Go to every component, just in case thery have to modify an initial condition (it is the case of the thermoviscoelastic elements)
	
	map<string,Body*,less<string> >::iterator itb;
	map<string,Joint*,less<string> >::iterator itj;
	
	/** Request the bodies */
	for(itb=mBodies.begin(); itb!=mBodies.end(); itb++)
	{
		itb->second->setInitialConditions(q,dq);
	}
	/** Request the joints */
	for(itj=mJoints.begin(); itj!=mJoints.end(); itj++)
	{
		itj->second->setInitialConditions(q,dq);
	}			
}

// ***

void MecSystem::update(const int& integratorType, const vector<double>& qConv, const vector<double>& dqConv)
{
	map<string,Body*,less<string> >::iterator i;
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) ((*j).second)->update(integratorType, qConv, dqConv);
	
	/** Request the bodies */
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		(*((*i).second)).update(integratorType, qConv, dqConv);
	}

	/** Request the joints */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		(*((*k).second)).update(integratorType, qConv, dqConv);
	}				
}

// ***

int MecSystem::setInertiaNumbers(int num, vector<int>& v)
{
	//Recursive	
	map<string,Point,less<string> >::iterator j;
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator kk;
	
    for(k=mChildSys.begin(); k!=mChildSys.end(); k++) num = ( ((*k).second)->setInertiaNumbers(num,v) );	/** Recursion */

	// Loop on its points
	for(j=mPoints.begin(); j!=mPoints.end(); j++)	num = ((*j).second).setInertiaNumbers(num,v);	

	// Loop on components (through bodies and joints)
	for(i=mBodies.begin(); i!=mBodies.end(); i++)	num = i->second->setInertiaNumbers(num,v);
	
	for(kk=mJoints.begin(); kk!=mJoints.end(); kk++)	num = kk->second->setInertiaNumbers(num,v);
	
	return(num);
}

// ***

void MecSystem::setTemperature(double& theta)
{
	map<string,Body*,less<string> >::iterator i;
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) j->second->setTemperature(theta);
	
	/** Request the bodies */
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		i->second->setTemperature(theta);
	}
	/** Request the joints */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		k->second->setTemperature(theta);
	}					
}

// ***

Body* MecSystem::getPTBodyLocal(string st)
{
	return mBodies[st];
}

// ***

void MecSystem::writeFileTemp(const double& time, ofstream & outF, string& bName)
{
	map<string,Body*,less<string> >::iterator i;
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) j->second->writeFileTemp(time,outF,bName);
	
	/** Request the bodies */
	/*
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		i->second->writeFileTemp(time,outF);
	}
	 */
	/*
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
			cout << "Nameeee = " << i->second->name << endl;
	}
*/
	mBodies[bName]->writeFileTemp(time,outF);
	/*
	for (int cont=1; cont<parentProc->vTempOut.size(); ++cont)
	{
		bName = parentProc->vTempOut[cont];
		cout << bName  << endl;
		mBodies[bName]->writeFileTemp(time,outF);
	}
	 */
	/*
	// Request the joints
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		k->second->writeFileTemp(time,outF);
	}
	*/
}

// ***

void MecSystem::setIniPosBody(string& iniJointName, vector<double>& iniValues)
{
	mBodies[iniJointName]->setIniPosBody(iniValues);
}

// ***

void MecSystem::setQatComponents(const vector<double>& q, const vector<double>& dq)
{
	map<string,Body*,less<string> >::iterator i;
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) j->second->setQatComponents(q,dq);
	
	/** Request the bodies */
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		i->second->setQatComponents(q,dq);
	}
	/** Request the joints */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		k->second->setQatComponents(q,dq);
	}						
}

// ***

void MecSystem::updateLagMul(void)
{
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) j->second->updateLagMul();
	
	/** Request the joints */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		k->second->updateLagMul();
	}							
}

// ***

void MecSystem::checkConvLagMul(bool& c)
{
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) j->second->checkConvLagMul(c);
	
	/** Request the joints */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		k->second->checkConvLagMul(c);
	}								
}

// ***

void MecSystem::calcProjMatrix(sparseMatrix& A, double& projPenal)
{
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) j->second->calcProjMatrix(A, projPenal);
	
	/** Request the joints */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		k->second->calcProjMatrix(A, projPenal);
	}									
}

// ***

void MecSystem::fillPhiq(vector<vector<double> >& Phiq, int& numd)
{
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) j->second->fillPhiq(Phiq,numd);
	
	/** Request the joints */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		k->second->fillPhiq(Phiq,numd);
	}											
}

// ***

void MecSystem::setSurfAtBodies()
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	
	Surface* surf;
	string surfName;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->setSurfAtBodies();	// Recursion
	
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		surfName = i->second->surfName;
		surf = parentProc->mSurfaces[surfName];
		i->second->setSurfAtBody(surf);
	}
}

// ***
/*
void MecSystem::setLocalPointsAtBodies()
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->setLocalPointsAtBodies();	// Recursion
	
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		if ( !(i->second->localPointsOption).compare("all") )
		{
			   i->second->setLocalPointsAtBody();
		}
	// Different mesh, TO BE DONE
	}
}
*/
// ***

/*
void MecSystem::setBodiesAtJoints()
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Joint*,less<string> >::iterator i;
	
	string bodyName1, bodyName2;
	Body* pBody1;
	Body* pBody2;
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->setBodiesAtJoints();	// Recursion
	
	for(i=mJoints.begin(); i!=mJoints.end(); i++)
	{
		bodyName1 = i->second->bodyName1;
		bodyName2 = i->second->bodyName2;
		if (bodyName1.compare("") && bodyName2.compare("") )
		{
			// Pointers to non-local bodies TO BE DONE
			pBody1 = getPTBodyLocal(bodyName1);
			pBody2 = getPTBodyLocal(bodyName2);
			i->second->setBodiesAtJoint(pBody1, pBody2);
		}
	}
}
*/

// ***

void MecSystem::createContacts()
{
	map<string,MecSystem*,less<string> >::iterator j;
	map<string,Joint*,less<string> >::iterator k;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) j->second->createContacts();
	
	/** Request the joints */
	for(k=mJoints.begin(); k!=mJoints.end(); k++)
	{
		k->second->createContacts();
	}											
}


// ***

void MecSystem::addNoDofsPoint(string& pname, string& stDofs, double& endtime )
{
	parentProc->addNoDofsPoint(pname,stDofs,endtime);
}

// ***

void MecSystem::repairChangedPointNames()
{
	map<string,MecSystem*,less<string> >::iterator j;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) j->second->repairChangedPointNames();
		
	string sysName;
	MecSystem* pSys;
	string oldName;
	string finalName;
	map<string,string,less<string> >::iterator it;
	map<string,Body*,less<string> >::iterator i;
	map<string,Joint*,less<string> >::iterator k;
	vector<double> coords;

	//cout << "En repairChangedNames de system " << name << " , has size " <<  mChangedPoints.size() <<  endl;

	for (it=mChangedPoints.begin(); it!=mChangedPoints.end(); ++it)
	{
		// Search through the systems for the final name
		oldName = it->first;
		//cout << "Searching for the new name of " << oldName << endl;
		sysName = getStringHead(it->second);
		pSys = getPTSystem(sysName);
		finalName = pSys->getNewName(it->second);
		//cout << "finalName is  " << finalName << endl;
		
		// Update mChangedPoints
		it->second = finalName;

		// Change names in system
		coords = pointNames[oldName];
		pointNames.erase(oldName);		
		pointNames[finalName] = coords;	
		
		// Change names in bodies and joints
		for(i=mBodies.begin(); i!=mBodies.end(); i++)
		{
			i->second->changePointName(oldName, finalName);
		}
		for(k=mJoints.begin(); k!=mJoints.end(); k++)
		{
			k->second->changePointName(oldName, finalName);
		}		
	}
}

// ***

string MecSystem::getNewName(string& oldName)
{
	map<string,string,less<string> >::iterator it;
	int num = mChangedPoints.count(oldName);
	if (num!=0)
	{	
		it = mChangedPoints.find(oldName);
		return it->second;
	}
	else return oldName;
}

// ***

void MecSystem::fillMSmooth(sparseMatrix& Ma, const vector<double>& q, double time)
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->fillMSmooth(Ma,q,time);	/** Recursion */
	
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		((*i).second)->fillMSmooth(Ma,q,time);
	}
}

// ***

void MecSystem::fillFSmooth(vector<vector<double> >& FSmooth, vector<double>& q_in, vector<double>& dq_in, double& time)
{
	map<string,MecSystem*,less<string> >::iterator k;
	map<string,Body*,less<string> >::iterator i;
	
	for(k=mChildSys.begin(); k!=mChildSys.end(); k++) ((*k).second)->fillFSmooth(FSmooth,q_in,dq_in,time);	/** Recursion */
	
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		((*i).second)->fillFSmooth(FSmooth,q_in,dq_in,time);
	}	
}

// ***

double MecSystem::getStressValue(int mm, int nn)
{
	if (mm<2)	return(parentProc->FSmooth[mm][nn]);
	else	return(parentProc->FSmoothDerived[mm-2][nn]);
}

// ***
				   
int MecSystem::checkComponents(int state)
{
	map<string,Body*,less<string> >::iterator i;
	map<string,MecSystem*,less<string> >::iterator j;
	
	/** Recursive */
	for(j=mChildSys.begin(); j!=mChildSys.end(); j++) state = ((*j).second)->checkComponents(state);
	
	/** Request the bodies */
	for(i=mBodies.begin(); i!=mBodies.end(); i++)
	{
		//state = (*((*i).second)).checkComponents(state);
		state = ((*i).second)->checkComponents(state);
	}
	return state;
}

// ***
  
void MecSystem::readTorsionFreeBeam(ifstream& in, string &auxs)
{
	string auxs1;
	string firstName = "";			/** Name of the point at the first end */
	string secondName = "";			/** Name of the point at the second end */
	vector<double> firstP;	/** Coords of first end point */
	vector<double> secondP;	/** Coords of second end point */
	string na = "";
	double ka, kd;
	string matName;
	double area;
	int nPieces;
	int beginNum;
	int i;
	int j;
	stringstream pointNum;
	string pointName;
	string ppointName;
	string midpointName;
	vector<double> coords(dim);
	int contNum = 0;
	vector<string> auxsv3(3);
	vector<vector<string> > auxsv;
	vector<string> allPointNames;
	double totalLength;
	double isec;
	
	// FIRST POINT
	in >> auxs1;	// P1 Label
	in >> firstName;		/** Read point name */
	firstP.resize(dim);
	if(checkForSlash(firstName)==false && !pointNames.count(firstName))	/** P1 belongs to this system because does not have slash and it is in the point list */
	{
	   addSystemName(firstName,name);
	}
	// Check if the point is local to the system or not
	if( checkIfPointSystem(firstName)==1 )
	{
	   firstP = pointNames[firstName];	// The point is local to the system
	}
	else												// The point is NOT local to the system
	{
	   na = getStringHead(firstName);	
	   MecSystem* pm = getPTSystem(na);
	   firstP = pm->pointNames[firstName];
	}
	
	// SECOND POINT
	in >> auxs1;	// P2 Label
	in >> secondName;		/** Read point name */
	secondP.resize(dim);
	if(checkForSlash(secondName)==false && !pointNames.count(secondName))	/** P2 belongs to this system because does not have slash and it is in the point list */
	{
		addSystemName(secondName,name);
	}
	// Check if the point is local to the system or not
	if( checkIfPointSystem(secondName)==1 )
	{
		secondP = pointNames[secondName];	// The point is local to the system
	}
	else											// The point is NOT local to the system
	{
		na = getStringHead(secondName);	
		MecSystem* pm = getPTSystem(na);
		secondP = pm->pointNames[secondName];
	}

	totalLength = lTwoNorm(firstP,secondP);
	
	in >> auxs1;				/** Reads the AREA label */
	in >> area;					/** Reads the area value */
	in >> auxs1;				// Reads the I label
	in >> isec;				 // Reads the I value
	in >> auxs1;				/** Reads the MAT label */
	in >> matName;				/** Reads the material name */
	in >> auxs1;				/** Reads the NPIECES label */
	in >> nPieces;				/**	Reads the npieces value */
	in >> auxs1;				/** Reads the BEGINNUM label */
	in >> beginNum;				/** Reads the beginnum value */
	//in >> auxs1;				// Reads the K label
	//in >> ka;					// Reads the ka value
	//in >> kd;					// Reads the kd value
	
	//cout << "nPieces = " << nPieces << endl;
	//cout << "firstName = " << firstName << endl;
	//cout << "secondName = " << secondName << endl;

	double factor = 1.0/(nPieces+1);
	
	allPointNames.push_back(firstName);
	
	if (nPieces==1)	// Just one piece
	{
		ppointName = firstName;
		pointName = secondName;
		pointNum << (beginNum + contNum);						/** Point number "int" is casted to "string" */
		pointNum >> midpointName;
		addSystemName(midpointName,name);	/** System name is added to the point name */
		for(j=0;j<dim;j++) coords[j] = (firstP[j]+secondP[j])/2.0;	/** Computes coordinates of point */
		addPointName(midpointName,coords);			/** Adds point name and coordinates at parent System */
		pointNum.clear();
//		contNum+=1;
		
		auxsv3[0] = ppointName;
		auxsv3[1] = midpointName;
		auxsv3[2] = pointName;
		auxsv.push_back(auxsv3);
		allPointNames.push_back(midpointName);
	}
	else // More than one piece
	{
		for (i=1; i<=nPieces; i++)
		{
			//cout << "i = " << i << endl;
			if (i==1)				// First
			{
				ppointName = firstName;
				
				pointNum << (beginNum + contNum);						/** Point number "int" is casted to "string" */
				pointNum >> midpointName;
				addSystemName(midpointName,name);	/** System name is added to the point name */
				for(j=0;j<dim;j++) coords[j] = (firstP[j]+factor*(secondP[j]-firstP[j]));	/** Computes coordinates of point */
				addPointName(midpointName,coords);							/** Adds point name and coordinates at parent System */
				pointNum.clear();
				contNum+=1;
				//
				pointNum << (beginNum + contNum);						/** Point number "int" is casted to "string" */
				pointNum >> pointName;
				addSystemName(pointName,name);	/** System name is added to the point name */
				for(j=0;j<dim;j++) coords[j] = (firstP[j]+2.*factor*(secondP[j]-firstP[j]));	/** Computes coordinates of point */
				addPointName(pointName,coords);							/** Adds point name and coordinates at parent System */
				pointNum.clear();
				contNum+=1;
				//
				allPointNames.push_back(midpointName);
				allPointNames.push_back(pointName);
			}
			else if (i==nPieces)	// Last
			{
				ppointName = midpointName;
				midpointName = pointName;
				
				//cout << "ppointName  = " << ppointName << endl;
				//cout << "midpointName  = " << midpointName << endl;
				pointName = secondName;
			}
			else					// Intermediate
			{
				ppointName = midpointName;
				midpointName = pointName;
				
				pointNum << (beginNum + contNum);						/** Point number "int" is casted to "string" */
				pointNum >> pointName;
				addSystemName(pointName,name);	/** System name is added to the point name */
				for(j=0;j<dim;j++) coords[j] = (firstP[j]+(i+1)*factor*(secondP[j]-firstP[j]));	/** Computes coordinates of point */
				addPointName(pointName,coords);							/** Adds point name and coordinates at parent System */
				pointNum.clear();
				contNum+=1;
				allPointNames.push_back(pointName);
			}
			auxsv3[0] = ppointName;
			auxsv3[1] = midpointName;
			auxsv3[2] = pointName;
			
			auxsv.push_back(auxsv3);
		}
	}
	allPointNames.push_back(secondName);
	
	mBodies[auxs] = new MassLinkChain(in,auxs,matName,area,allPointNames,this);

	mJoints[auxs] = new TorsionFreeBeam(in,auxs+"c",this, allPointNames, auxsv, matName, totalLength, area, isec);
}

// ***

void MecSystem::readRigidBar(ifstream& in, string &auxs)
{
	string auxs1;
	vector<string> nms(2);
	double masa;
	
	// FIRST POINT
	in >> auxs1;	// P1 Label
	in >> nms[0];		/** Read point name */
	if(checkForSlash(nms[0])==false && !pointNames.count(nms[0]))	/** P1 belongs to this system because does not have slash and it is in the point list */
	{
		addSystemName(nms[0],name);
	}
	
	// SECOND POINT
	in >> auxs1;	// P2 Label
	in >> nms[1];		/** Read point name */
	if(checkForSlash(nms[1])==false && !pointNames.count(nms[1]))	/** P2 belongs to this system because does not have slash and it is in the point list */
	{
		addSystemName(nms[1],name);
	}
	// MASS
	in >> auxs1;	// MASS Label
	in >> masa;		/** value */
	
	mBodies[auxs] = new MassLinkChain(in,auxs,masa,nms,this);
	mJoints[auxs] = new RigidLink(in,auxs,nms,this);
}
